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1. Introduction 

At the Geophysical Fluid Dynamics Laboratory, research is carried out for 
the purpose of understanding various aspects of climate, such as its 
variability, predictability, stability and sensitivity. The atmosphere and 
oceans are modelled mathematically and their phenomenology studied by computer 
simulation methods. The present paper will discuss the present state-of-the- 
art in the computer simulation of large scale oceans on the CYBER 205. While 
atmospheric modelling differs in some aspects, the basic approach used is 
similar. 

The equations of the ocean model will be presented in the following 
section along with a short description of the numerical techniques used to find 
their solution. Section 3 will deal with computational considerations and a 
typical solution will be presented in section 4- 


2. Equations of the model 

The model presented here is the multilevel numerical model described in 
Bryan (1969). The continuous equations will be given. A detailed description 
of the finite difference formulation may be found in the above work. The 
equations of motion are the Navier-Stokes equations written in spherical 
coordinates and modified by the Boussinesq approximation. Let m=sec0, 
n=sin^, u=aXm ^ and v=a0, where a is the radius of the earth, 

0 the latitude and X the longitude. It is convenient to define the 
advect ion operator 

r(^>=»a"l [(u^)x+(v^»"l)0] +(w^)2. (1) 

The equations of motion on a sphere are 

ut*r<u)-2nnv=-M‘^(P/Po^)x**^^' 
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Vt*r<v)*2nnu= 


<3) 


r(i)=o. 


(4) 


gp=-Pz' 


(5) 


where Pq is unity in cgs units. The conservation equations for the 
temperature and salinity are 

Tt+r<T)=F^ (g 

Stores) (7 

The terms in F contain effects of nixing as well as external driving forces. 
The equation of state 


p=P(T,S,z) (8) 

is an empirically derived formula relating the local density of seawater to 
temperature^ salinity and depth. 

The set of equations (1-8) are cast into finite difference form. The 
prognostic equations (2, 3,6,7) are solved as an initial value problem, placing 
all terms except the local time derivative on the right hand side and carrying 
out timesteps to predict new values of velocity, temperature and salinity on a 
prescribed mesh covering the model ocean domain. Given a certain configuration 
of steady wind driving and differential surface heating (both entering through 
the F terms), a statistical steady state is approached asymptotically in time. 
Time scale analysis of Eqs.(6,7) reveals that 0(1000) years of integration is 
needed to bring the sluggish abyssal layers of the ocean model into a steady 
state . 


28 



Ill 


3. Computational considerations 

Let us consider a rectangular ocean basin model comparable in size to the 

N- Atlantic Ocean. It extends 60® in longitude, 65® in latitude and 

4000 meters in depth. It is desirable to cover this domain with a mesh fine 

enough to resolve mesoscale <0(100 km)) eddies which play an important role in 

transporting various properties through the ocean. The minimum resolution 
needed for this purpose is roughly l/3rd degree in latitude and somewhat 
larger, say .4 degree in longitude due to the convergence of meridians on the 
globe. This results in a horizontal grid space of 150x195 points. Vertically, 
18 levels are needed to resolve the scales of interest. This brings the total 
to just over 1/2 million grid points for which Eqs.(l~8) must be evaluated each 
timestep. 

The longest timestep which can be used without incurring numerical 
instability is given by the Courant*Fr iedrichs-Lewy condition 

c^t/^<l (9) 


where c is the phase velocity of the fastest moving wave in the ocean. Since 
high speed external gravity waves have been filtered from this model by the 
condition w-0 at the surface, the fastest wave is that associated with the 
internal density gradients (internal gravity wave) which has a speed of roughly 
3m/sec. The smallest Ax occurs at the northern wall of the model due to 
convergence of meridians, and is about 20 km. The resultant At is such 
that roughly 5000 timesteps are necessary to integrate one year. Therefore, 5 
million timesteps, or 2. 5x10^2 grid point evaluations of Eqs.(l-8), are 
required to integrate this model to a steady state- Even the fastest modern 
day computers cannot accomplish this task in a reasonable time, although steady 
progress is being made. The former computer at GFDL, the Texas Instruments 
ASC, took 15 seconds to compute one time step on the above model. At this 
speed, 2.4 years of computing would be needed to reach a steady state solution. 
Clearly, compromises must be made in designing experiments which are achievable 
in a reasonable amount of computer time. This may involve reducing the domain 
size, or integrating for a shorter period, or both. (Interesting results may be 
obtained from an integration of 0(10) years, particularly for the upper ocean 
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where ti»e scales of adjustnent are relatively short.) The greater the 
computational speed which can be attained, the less severe the compromises must 
be. 

In converting the ASC ocean model to the CYBER 205, the most fundamental 
alteration of the code had to do with the treatment of land masses- 
Previoualy, the computation was carried out only over ocean points by making 
the DO loop limits functions of the placement of land. The contiguity 
requirement of the 205 for vectorization allows only the innermost of the three 
dimensional loops to vectorize in this case. An alternative method of handling 
land is to compute all points as if they were ocean and, at the end of the 
timestep, restore the land to its specified value using a masking array. 
Contiguity is then satisfied and vectorization is enabled through two 
dimensions. (The third dimension cannot be vectorized because it is cycled 
through memory from disc.) By using the latter technique, the typical vector 
length in the computation is increased from 150 in the example above (east-west 
dimension) to 2700 (east-west times depth dimension) resulting in a 
considerable decrease in the relative time spent in vector startup. 

An additional time saving has been accomplished in an area of the code 
which is used heavily, but is inherently unvectorizable due to a recursive 
property. Using 08 calls to insert machine language directly into the FORTRAN, 
CDC personnel have “unrolled" this loop, greatly improving on the code 
generated by the compiler for the equivalent FORTRAN loop. 

The use of half -precision on all floating point variables has resulted in 
a gain of only about 15X in overall running speed, although sections of the 
code which are 100^ vectorized increase in speed by roughly 405i. Additional 
work is needed to determine why the overall gain is so small considering the 
high degree of vectorization of the code. 

Since the model above is too large to fit into core memory entirely, data 
is cycled through memory from disc as it is needed each timestep. If this disc 
transfer cannot be buffered sufficiently well, computation ceases while waiting 
for the I/O to finish. The result is that the computer may not be used 
efficiently, particularly if the other jobs running concurrently have the same 
difficulty. Until recently, this was a severe problem on the 205. The above 
model, when in the 205 alone, ran only about 15^ of the wall clock time. 

Improved I/O schemes have been developed by CDC personnel at GFDL and currently 
the same model runs about SOX of the wall clock time when alone. This compares 
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favorably with I/O efficiencies on the ASC- 

The CYBER 205 version of the model described above currently takes 4 
seconds to compute one timestep, almost a factor of 4 faster than the ASC. 

While this speed still does not make the experiment proposed at the beginning 
of this section feasible, the compromises which are necessary to produce an 
attainable solution are much less severe than before- One such experiment will 
be described in the following section. 

4. An ocean simulation experiment 

If one wishes to study the effects of topography on the dynamics of the 
Gulf Stream, an argument can be made that it is not necessary to consider a 
domain as large as the one proposed earlier, and that several decades of 
integration is sufficient. Therefore, let us reduce the domain from 65 to 27 
degrees in latitude and from 60 to 32 degrees in longitude- Also, for this 
purpose, the vertical resolution may be decreased from 18 layers to 5 layers. 
This produces a model which takes approximately one hour of 205 time to inte- 
grate one year of ocean time. Applying surface wind stress and differential 
heating similar to that of the N. Atlantic, this model has been integrated from 
rest a total of 20 years. The resulting temperature pattern at the second 
layer, centered at 212 meters depth, is shown in Fig. 1. The land mass in the 
northwest corner simulates the gross features of the U-S. east coast. A conti- 
nental shelf and slope is also included in this solution. The simulated Gulf 
Stream is revealed by the tightly packed isotherms along the coast and bending 
out to sea at the point representing Cape Hatteras. In agreement with 
observations, there exist both cold and warm core "rings'* which have broken 
from the Stream and are drifting westward. An example of the former is 
centered at about 70«W, 30«N and of the latter at GS^W, 37«N. 

Three other experiments have been carried out in this series, altering the 
topography along the western boundary to study its effect on the path and 
behavior of the Gulf Stream. 
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